## Agroforestry Paper Analysis
## Author : l.orero@cgiar.org
## This is an script to analyze and create tables and graphs for paper
## Last update: 27th May 2021

# Create working datasets import
# tree dataset
library(readxl)
nyando1 <- read_excel("C:/Users/LORERO/OneDrive - CGIAR/Desktop/ICRAF/Current Projects/Nyando Project/Data/2019 Nyando Data.xlsx", 
    sheet = "set")
View(nyando1)

nyando1$Block = as.factor(nyando1$Block)
nyando1$type_r = as.factor(nyando1$type_r)

# Create subset for Total Lower
# Project
nyandolp <- subset(nyando1, type_r == 
    "Lower Project")
View(nyandolp)

# Create subset for Total Lower
# Control
nyandolc <- subset(nyando1, type_r == 
    "Lower Control")
View(nyandolc)

# Create subset for Total Middle
# Project
nyandomp <- subset(nyando1, type_r == 
    "Middle Project")
View(nyandomp)

# Create subset for Total Middle
# Control
nyandomc <- subset(nyando1, type_r == 
    "Middle Control")
View(nyandomc)

# Create a population DB for all
# data
library(dplyr)
library(tidyr)
library(tidyverse)

steph <- nyando1
steph$count <- 1
steph1 <- steph[, c("IntervieweeName", 
    "Block", "type_r", "ha", "count")]
steph2 <- aggregate(. ~ IntervieweeName + 
    Block + type_r, steph1, sum)

# Create 5 DBH Classes of five DBH
# classes of [<10, 10-20, 20-30,
# 30-40, and >40 cm]

## Tree density of less than 10 cm
## (<10cm)###########################
nyando1_d1 <- filter(nyando1, DBH < 
    10)
nyando1_d1$count <- 1
nyando_d11 <- nyando1_d1[, c("IntervieweeName", 
    "Block", "type_r", "ha", "count")]
aggdata1 <- aggregate(. ~ IntervieweeName + 
    Block + type_r, nyando_d11, sum)
aggdata1$ha <- steph2$ha[match(aggdata1$IntervieweeName, 
    steph2$IntervieweeName)]
aggdata1$tree_d <- aggdata1$count/aggdata1$ha

## Lower vs Middle ANOVA
aov.1 <- aov(tree_d ~ Block, data = aggdata1)
summary(aov.1)

oneway.test(tree_d ~ Block, data = aggdata1)  #ANOVA Unequal var
## One-way analysis of means (not
## assuming equal variances)

oneway.test(tree_d ~ type_r, data = subset(aggdata1, 
    Block == "Middle"))
## One-way analysis of means (not
## assuming equal variances)

oneway.test(tree_d ~ type_r, data = subset(aggdata1, 
    Block == "Lower"))
## One-way analysis of means (not
## assuming equal variances)

# 10-20cm
nyando_d2 <- filter(nyando1, DBH >= 
    10 & DBH <= 20)
nyando_d2$count <- 1
nyando_d21 <- nyando_d2[, c("IntervieweeName", 
    "Block", "type_r", "hectare", "count")]
aggdata2 <- aggregate(. ~ IntervieweeName + 
    Block + type_r, nyando_d21, sum)
aggdata2$ha <- steph2$hectare[match(aggdata2$IntervieweeName, 
    steph2$IntervieweeName)]
aggdata2$tree_d <- aggdata2$count/aggdata2$ha

## Lower vs Middle ANOVA
aov.2 <- aov(tree_d ~ Block, data = aggdata2)
summary(aov.2)

oneway.test(tree_d ~ Block, data = aggdata2)  #ANOVA Unequal var

oneway.test(tree_d ~ type_r, data = subset(aggdata2, 
    Block == "Middle"))

oneway.test(tree_d ~ type_r, data = subset(aggdata2, 
    Block == "Lower"))

# 20-30CM
nyando_d3 <- filter(nyando1, DBH > 20 & 
    DBH <= 30)
nyando_d3$count <- 1
nyando_d31 <- nyando_d3[, c("IntervieweeName", 
    "Block", "type_r", "hectare", "count")]

aggdata3 <- aggregate(. ~ IntervieweeName + 
    Block + type_r, nyando_d31, sum)
aggdata3$ha <- steph2$hectare[match(aggdata3$IntervieweeName, 
    steph2$IntervieweeName)]

aggdata3$tree_d <- aggdata3$count/aggdata3$ha


## Lower vs Middle ANOVA
aov.3 <- aov(tree_d ~ Block, data = aggdata3)
summary(aov.3)

oneway.test(tree_d ~ Block, data = aggdata3)  #ANOVA Unequal var

## Control vs Project ANOVA
aov.3b <- aov(tree_d ~ type_r, data = aggdata3)
summary(aov.3b)

# TukeyHSD
TukeyHSD(aov.3b)

## Control vs Project ANOVA (Relaxing
## the homogeneity of variance
## assumption)
oneway.test(tree_d ~ type_r, data = subset(aggdata3, 
    Block == "Middle"))

oneway.test(tree_d ~ type_r, data = subset(aggdata3, 
    Block == "Lower"))

# 30-40cm
nyando_d4 <- filter(nyando1, DBH > 30 & 
    DBH <= 40)
nyando_d4$count <- 1
nyando_d41 <- nyando_d4[, c("IntervieweeName", 
    "Block", "type_r", "hectare", "count")]

aggdata4 <- aggregate(. ~ IntervieweeName + 
    Block + type_r, nyando_d41, sum)
aggdata4$ha <- steph2$hectare[match(aggdata4$IntervieweeName, 
    steph2$IntervieweeName)]

aggdata4$tree_d <- aggdata4$count/aggdata4$ha

## Lower vs Middle ANOVA
aov.4 <- aov(tree_d ~ Block, data = aggdata4)
summary(aov.4)

oneway.test(tree_d ~ Block, data = aggdata4)  #ANOVA Unequal var
## 

## Control vs Project ANOVA
aov.4b <- aov(tree_d ~ type_r, data = aggdata4)
summary(aov.4b)

# TukeyHSD
TukeyHSD(aov.4b)

## Control vs Project ANOVA (Relaxing
## the homogeneity of variance
## assumption)
oneway.test(tree_d ~ type_r, data = subset(aggdata4, 
    Block == "Middle"))

oneway.test(tree_d ~ type_r, data = subset(aggdata4, 
    Block == "Lower"))
## 

# >40cm
nyando_d5 <- filter(nyando1, DBH > 40)
nyando_d5$count <- 1
nyando_d51 <- nyando_d5[, c("IntervieweeName", 
    "Block", "type_r", "hectare", "count")]

aggdata5 <- aggregate(. ~ IntervieweeName + 
    Block + type_r, nyando_d51, sum)
aggdata5$ha <- steph2$hectare[match(aggdata5$IntervieweeName, 
    steph2$IntervieweeName)]

aggdata5$tree_d <- aggdata5$count/aggdata5$ha

## Lower vs Middle ANOVA
aov.5 <- aov(tree_d ~ Block, data = aggdata5)
summary(aov.5)

oneway.test(tree_d ~ Block, data = aggdata5)  #ANOVA Unequal var
## 

## Control vs Project ANOVA
aov.5b <- aov(tree_d ~ type_r, data = aggdata5)
summary(aov.5b)

# TukeyHSD
TukeyHSD(aov.5b)

## Control vs Project ANOVA (Relaxing
## the homogeneity of variance
## assumption)
oneway.test(tree_d ~ type_r, data = subset(aggdata5, 
    Block == "Middle"))
## 

# One way ANOVA DBH1 Lower vs Middle
db1 <- aov(DBH ~ type_r, data = nyando1_d1)
summary(db1)

nyando_c <- nyando1
nyando_c$count <- 1
nyando_c1 <- nyando_c[, c("IntervieweeName", 
    "Block", "type_r", "hectare", "count")]

aggdatac <- aggregate(. ~ IntervieweeName + 
    Block + type_r, nyando_c1, sum)
aggdatac$ha <- steph2$hectare[match(aggdatac$IntervieweeName, 
    steph2$IntervieweeName)]
aggdatac$tree_d <- aggdatac$count/aggdatac$ha

## Lower vs Middle ANOVA
aov.c <- aov(tree_d ~ Block, data = aggdatac)
summary(aov.c)

oneway.test(tree_d ~ Block, data = aggdatac)  #ANOVA Unequal var


## Control vs Project ANOVA
aov.cb <- aov(tree_d ~ type_r, data = aggdatac)
summary(aov.cb)

# TukeyHSD
TukeyHSD(aov.cb)

# Levene's Test for Homogeneity of
# Variance (center = median)
oneway.test(tree_d ~ type_r, data = aggdatac)

# AGB
# Calculations############################################
# AGB Calculations Combined sample
# of all DBH
nyando_gb <- nyando1
nyando_gb$num <- rep(1:length(nyando1$IntervieweeName))
nyando_agb <- nyando_gb[, c("Block", 
    "type_r", "Agroforestry Practice", 
    "AGB (Mg)", "num")]
aggdata_gb <- aggregate(. ~ num + Block + 
    type_r + `Agroforestry Practice`, 
    nyando_agb, sum)
aggdata_gb$agb = aggdata_gb$`AGB (Mg)`
aggdata_gb$agfp = aggdata_gb$`Agroforestry Practice`
aggdata_gb$Block = as.factor(aggdata_gb$Block)
aggdata_gb$type_r = as.factor(aggdata_gb$type_r)

## Lower vs Middle ANOVA for AGB
oneway.test(agb ~ Block, data = aggdata_gb)  #ANOVA Unequal var

## Lower Project AGB
aggdata_gb1 <- filter(aggdata_gb, type_r == 
    "Lower Project")

# Total AGB under Boundary Planting
sum(filter(aggdata_gb1, agfp == "Boundary planting")$agb)


# Total AGB under Riparian buffers
sum(filter(aggdata_gb1, agfp == "Riparian buffers")$agb)
## [1] 0.09491625

## Lower Control AGB
aggdata_gb2 <- filter(aggdata_gb, type_r == 
    "Lower Control")

## Middle Project AGB
aggdata_gb3 <- filter(aggdata_gb, type_r == 
    "Middle Project")

# Middle Control AGB
aggdata_gb4 <- filter(aggdata_gb, type_r == 
    "Middle Control")
